* Reset settings and initialize log file
launch, path("share/edexits")

*-------------------------------------------------------------------------------
* Price and Wasserman (2024), "The Summer Drop in Female Employment"
*
* Description: Show gender differences in hazard rates of transitioning from
*              the education sector to non-employment.
*-------------------------------------------------------------------------------


* Run specifications
*-------------------------------------------------------------------------------

if "$estimate" != "0" {
	* Load data on adult individuals
	gzuse "$basepath/data/derived/cps_bms_sample.dta.gz", clear
	keep pid tm month wtraked tmspline* weeks female empstat emp school occ1990

	* Record lagged employment status and occupation
	gen byte L_empstat = L.empstat
	gen int L_occ1990 = L.occ1990

	* Restrict to individuals employed in education last month
	keep if L.emp == 1 & L.school == 1 & !missing(wtraked)

	* Construct indicators for select transitions, conditional on being at risk
	gen byte haz_nonemp = (emp == 0)
	gen byte haz_nohours = (empstat != 10) if L_empstat == 10

	* Tag select occupations
	gen byte occ_code = .
	replace occ_code = 1 if L_occ1990 == 156
	replace occ_code = 2 if L_occ1990 == 157
	replace occ_code = 3 if L_occ1990 == 014
	replace occ_code = 4 if L_occ1990 == 808

	label define occ_code_lbl 0 "All education workers", replace
	label define occ_code_lbl 1 "Primary school teachers", add
	label define occ_code_lbl 2 "Secondary school teachers", add
	label define occ_code_lbl 3 "Managers in education", add
	label define occ_code_lbl 4 "School bus drivers", add
	label values occ_code occ_code_lbl

	* Create a pooled category as well
	expand = 2
	bysort pid tm: replace occ_code = 0 if _n == 2
	drop if missing(occ_code)

	* Aggregate to the sex x job-type level
	gcollapse (mean) haz_nonemp haz_nohours (first) month weeks tmspline* (rawsum) wtraked [pw = wtraked], by(female occ_code tm)

	* Store the sample
	tempfile core
	save `core'

	* Loop over sex
	foreach f of numlist 0 1 {
		* Reload the sample
		use if female == `f' using `core', clear
		tsset occ_code tm

		* Run a specification for each outcome x job type
		foreach v in "nonemp" "nohours" {
			foreach j of numlist 0/4 {
				quietly ivreg2 haz_`v' ib5.month D.tmspline* D.weeks if occ_code == `j' [aw = wtraked], bw($bandwidth) robust small
				process_estimates, path("edexits") model("f`f'_`v'_j`j'")
			}
		}
	}
}


* Prepare estimates
*-------------------------------------------------------------------------------

* Load estimates into memory
load_estimates, path("edexits")

* Prepare coefficient labels for each month
make_coeflabels


* Plot separation hazards from education jobs (left) and specific occs (right)
*-------------------------------------------------------------------------------

* Loop over outcome
foreach v in "nonemp" "nohours" {
	* Specify grayscale or color
	if "`v'" == "nonemp" {
		local col1 black
		local col2 gs6
		local mfcolor2 mfcolor(white)
	}
	else if "`v'" == "nohours" {
		local col1 "$col1"
		local col2 "$col2"
		local mfcolor2
	}

	* Left panel
	if "`v'" == "nonemp" {
		local ymin = -1.1
		local ymax = 6.0
		local ylabel "-1(1)6"
	}
	else if "`v'" == "nohours" {
		local ymin = -5
		local ymax = 40
		local ylabel "-5(5)40"
	}

	clear
	set obs 14
	gen rmin = `ymin'
	gen rmax = `ymax'
	gen rval = _n - 0.5

	#delimit ;
	coefplot
		(f1_`v'_j0, offset(-.05) recast(connected) color(`col1') ciopts(color(`col1')) label("Women"))
		(f0_`v'_j0, offset(+.05) recast(connected) color(`col2') ciopts(color(`col2')) `mfcolor2' label("Men")),
		coeflabels(`coeflabels')
		title("{bf:All workers in the education sector}")
		xtitle("")
		xlabel(, alternate)
		ytitle("Difference relative to May (p.p.)")
		yscale(range(`ymin' `ymax'))
		ylabel(`ylabel')
		graphregion(margin(small))
		plotregion(margin(vsmall))
		addplot(
			rarea rmax rmin rval if inrange(rval, 0.5, 5.5), color($ltgs) plotregion(margin(t=0 b=0)) below ||
			scatteri 0 0.5 0 13.5, recast(line) color(black) below)
		vertical
		legend(rows(1) size(*0.8))
		rescale(100);
	#delimit cr

	tempfile left
	graph save "`left'"

	* Right panel
	if "`v'" == "nonemp" {
		local ymin = -3
		local ymax = 15.5
		local ylabel "-3(3)15"
	}
	else if "`v'" == "nohours" {
		local ymin = -10
		local ymax = 60
		local ylabel "-10(10)60"
	}

	clear
	set obs 14
	gen rmin = `ymin'
	gen rmax = `ymax'
	gen rval = _n - 0.5

	foreach j of numlist 1/4 {
		if `j' == 1	local jlbl "Primary school teachers"
		if `j' == 2	local jlbl "Secondary school teachers"
		if `j' == 3	local jlbl "Managers in education"
		if `j' == 4	local jlbl "School bus drivers"

		#delimit ;
		coefplot
			(f1_`v'_j`j', keep(coef0 coef1 coef2 coef3 coef4) offset(-.05) recast(connected) color(`col1') ciopts(color(`col1')))
			(f0_`v'_j`j', keep(coef0 coef1 coef2 coef3 coef4) offset(+.05) recast(connected) color(`col2') ciopts(color(`col2')) `mfcolor2'),
			coeflabels(coef0 = "May" coef1 = "J" coef2 = "J" coef3 = "A" coef4 = "Sept.")
			title("{bf:`jlbl'}")
			xtitle("")
			xscale(range(0.5 5.5))
			ytitle("")
			yscale(range(`ymin' `ymax'))
			ylabel(`ylabel')
			yline(0)
			addplot(
				rarea rmax rmin rval if inrange(rval, 0.5, 5.5), color($ltgs) plotregion(margin(t=0 b=0)) below ||
				scatteri 0 0.5 0 5.5, recast(line) color(black) below)
			vertical
			legend(off)
			rescale(100);
		#delimit cr

		tempfile g`j'
		graph save "`g`j''"
	}

	* Create pdf file
	grc1leg "`g1'" "`g2'" "`g3'" "`g4'", rows(2) graphregion(margin(small)) plotregion(margin(vsmall)) imargin(vsmall)
	tempfile right
	graph save "`right'"

	grc1leg "`left'" "`right'", rows(1) imargin(vsmall)
	nicepdf "$basepath/output/edexits_`v'.pdf", panels(2) cropbuffer(1) indirect replace
}

* Close the log file
unlaunch
